Impianto di valvole cardiache: regressioni multiple di Cox a confronto

Sommario
Lo studio osservazionale che ha permesso la raccolta dei dati heart_valve riguarda pazienti che hanno subito un impianto di valvole cardiache. Si vuole dunque indagare e modellare la sopravvivenza di questi pazienti, tramite una raccolta dati durata oltre 10 anni.
In particolare, si esamina il ruolo dell’indice di massa ventricolare sinistra sulla prognosi. A questo scopo si sono costruiti due modelli di regressione multipla di Cox, uno base e uno aumentato. Il modello aumentato differisce da quello base per il solo inserimento della variabile di indice di massa ventricolare sinistra.
Successivamente alla validazione dei modelli e allo studio delle loro assunzioni di base, si è studiata la loro capacità predittiva del rischio, con un’attenzione particolare al beneficio che un modello aumentato può portare.

1 Dataset e preprocessing

Variabili incluse nel dataset:

• Paz.id: identificativo del paziente (stringa numerica)

• log.lvmi: logaritmo naturale dell’indice di massa ventricolare sinistra misurato al basale. Questa variabile si presenta standardizzata.

• fuyrs: tempo di follow-up dalla chirurgia (anni).

• status: indicatore di evento (1 = morto; 0 = perso al follow up).

• sex: genere del paziente (0 = M; 1 = F).

• age: età del paziente alla chirurgia (anni).

• con.cabg: presenza concomitante di bypass coronarico (1 = si; 0 = no).

• creat: creatinina serica pre-operatoria (\(\frac{\mu \text{mol}}{\text{mL}}\)).

• lv: frazione di eiezione ventricolare sinistra pre-operatoria (1 = buona, 2 = moderata, 3 = scarsa).

• sten.reg.mix: emodinamica della valvola aortica (1 = stenosi, 2 = rigurgito, 3 = misto).

Cominciamo importando i pacchetti.

# Stimatore di Kaplan-Meier
library(prodlim)

# Stima dei modelli di Cox
library(survival)

# Librerie grafiche
library(ggplot2)
library(corrplot)

# Inserimento di splines nei modelli
library(splines)

# Validazione dei modelli
library(riskRegression)

# Net Benefit Plot
library(dcurves)

# funzione plotHR
library(Greg)

# Clusters
library(cluster)

Dopo avere importato i dati, esaminiamo i primi valori riportati nelle colonne.

heart_valve <- read.table("data.txt",
     na.strings = ".",
     header = T, row.names = NULL
)
str(heart_valve)

## 'data.frame':    255 obs. of  10 variables:
##  $ paz.id      : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ log.lvmi    : num  4.78 4.74 5.26 4.23 5.14 ...
##  $ fuyrs       : num  4.96 9.66 7.92 4.04 8.82 ...
##  $ status      : int  0 0 0 0 0 1 0 0 0 0 ...
##  $ sex         : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ age         : num  75.1 45.8 63.3 61.4 53.6 ...
##  $ con.cabg    : int  1 0 0 1 0 1 0 0 0 0 ...
##  $ creat       : int  103 76 130 92 109 111 74 136 76 87 ...
##  $ lv          : int  1 2 1 2 1 1 2 1 1 2 ...
##  $ sten.reg.mix: int  1 1 2 2 2 1 3 2 1 3 ...

Notiamo che alcune variabili non sono state automaticamente inserite come discrete, quindi proseguiamo manualmente.

heart_valve$paz.id <- as.character(heart_valve$paz.id)
heart_valve$sex <- factor(heart_valve$sex)
heart_valve$con.cabg <- factor(heart_valve$con.cabg)
heart_valve$lv <- factor(heart_valve$lv)
heart_valve$sten.reg.mix <- factor(heart_valve$sten.reg.mix)

Inoltre vengono ricodificate in livelli le variabili continue. Questo ci permetterà di confrontarle con variabili discrete. Verranno scartate a favore di quelle originali, quando invece faranno parte dei regressori dei modelli.

heart_valve$age_lev <- cut(heart_valve$age, breaks = c(0, 50, 65, 100))
heart_valve$creat_lev <- cut(heart_valve$creat, breaks = c(0, 80, 120, 300))
heart_valve$log.lvmi_lev <- cut(heart_valve$log.lvmi, breaks = c(3.5, 4.8, 5.2, 6.5))

2 Analisi descrittive

Proseguiamo ora con le analisi descrittive dei dati, partendo da semplici indici statistici riguardo le distribuzioni, e continuando con la loro visualizzazione.

summary(heart_valve)

##     paz.id             log.lvmi         fuyrs              status       sex    
##  Length:255         Min.   :3.556   Min.   : 0.02192   Min.   :0.0000   0:182  
##  Class :character   1st Qu.:4.763   1st Qu.: 3.38904   1st Qu.:0.0000   1: 73  
##  Mode  :character   Median :5.043   Median : 5.03288   Median :0.0000          
##                     Mean   :5.055   Mean   : 5.32429   Mean   :0.2118          
##                     3rd Qu.:5.289   3rd Qu.: 7.38219   3rd Qu.:0.0000          
##                     Max.   :6.352   Max.   :10.73973   Max.   :1.0000          
##       age        con.cabg     creat       lv      sten.reg.mix     age_lev   
##  Min.   :23.33   0:177    Min.   : 50.0   1:146   1:158        (0,50]  : 27  
##  1st Qu.:58.11   1: 78    1st Qu.: 81.0   2: 87   2: 50        (50,65] : 76  
##  Median :67.84            Median : 96.0   3: 22   3: 47        (65,100]:152  
##  Mean   :65.94            Mean   :102.1                                      
##  3rd Qu.:75.14            3rd Qu.:116.5                                      
##  Max.   :89.12            Max.   :264.0                                      
##      creat_lev      log.lvmi_lev
##  (0,80]   : 61   (3.5,4.8]:73   
##  (80,120] :140   (4.8,5.2]:97   
##  (120,300]: 54   (5.2,6.5]:85   
##                                 
##                                 
##

Notiamo che quasi l’\(80\%\) dei soggetti in analisi sono persi al follow up, mentre il restante \(20\%\) circa presenta l’evento: non sono quindi presenti soggetti censurati. Tenendo in considerazione il lungo periodo durante il quale i pazienti vengono seguiti e il fatto che soggetti che non presentano problemi post-operatori non vengano probabilmente controllati di frequente, nel modellare la probabilità di sopravvivenza si assumono i soggetti persi al follow-up come right censored.

ggplot(
     heart_valve,
     aes(x = as.factor(status), fill = as.factor(status))
) +
     geom_bar() +
     scale_fill_hue(c = 40) +
     theme(legend.position = "none") +
     ggtitle("Distribuzione di evento") +
     xlab("Evento") +
     ylab("Frequenza assoluta")

Vediamo che il tempo di follow-up dei pazienti ha una distribuzione approssimativamente normale attorno alla media di 5 anni.

hist(heart_valve$fuyrs,
     xlab = "Anni", ylab = "Frequenza assoluta",
     main = "Distribuzione del tempo di follow-up"
)

Si osservi che più della metà dei pazienti presentano una valvola aortica stenotica e una frazione di eiezione ventricolare sinistra nella norma. Vi è anche una prevalenza di uomini nel campione. Questo è spiegabile dal fatto che, supponendo di avere un campione significativo, i problemi cardiaci e i conseguenti interventi alla valvola sono molto più predominanti negli uomini rispetto alle donne. Inoltre, circa 2/3 dei pazienti ha ricevuto un intervento con presenza concomitante di bypass coronario.

par(mfrow = c(2, 2))

barplot(table(heart_valve$sex),
     main = "Genere paziente"
)
barplot(table(heart_valve$con.cabg),
     main = "Presenza bypass coronarico"
)
barplot(table(heart_valve$lv),
     main = "Fraz. eiezione ventricolare sinistra"
)
barplot(table(heart_valve$sten.reg.mix),
     main = "Emodinamica valvola aortica"
)

Come atteso l’età dei pazienti si concentra sulle fasce di popolazione più anziane, trattandosi di interventi cardiaci. Il logaritmo del lvmi ha una distribuzione approssimativamente normale attorno alla media con alcuni valori estremi inferiori a 4 o superiori a 6. Il livello di creatinina serica pre-operatoria ha invece una distribuzione asimmetrica con una coda a destra. I valori standard di creatina serica vanno tra i 65 e i 130 microgrammi al millilitro, pertanto si riconoscono numerosi valori ben oltre le soglie adeguate, con individui che superato addirittura i 200 microgrammi al millilitro.

par(mfrow = c(1, 3))

hist(heart_valve$log.lvmi,
     xlab = "Left Ventricular Mass Index",
     main = "Log Left Ventricular Mass Index"
)
hist(heart_valve$age,
     xlab = "Anni",
     main = "Età del paziente"
)
hist(heart_valve$creat,
     xlab = "Creatina serica",
     main = "Creatinina serica"
)

Ora, sempre a livello descrittivo, viene riprodotta la funzione di sopravvivenza nel tempo per la totalità dei soggetti. Si evince che la sopravvivenza di un paziente che ha subito un impianto di valvole cardiache a 5 anni dall’operazione è dell’\(80\%\), mentre dopo 10 anni è poco superiore al \(50\%\).

fit.surv <- prodlim(Hist(fuyrs, status) ~ 1,
     data = heart_valve
)
plot(fit.surv)

Studiando ora la presenza concomitante di bypass coronario, scopriamo che questo porta ad una sopravvivenza significativamente inferiore sulla quasi totalità del periodo di 10 anni.

fit.surv_cc <- prodlim(Hist(fuyrs, status) ~ con.cabg,
     data = heart_valve
)
plot(fit.surv_cc, legend = FALSE)

Anche la frazione di eiezione ventricolare sinistra pre-operatoria sembra essere un fattore di forte impatto sulla sopravvivenza degli individui successivamente all’inserimento della alla valvola. Risulta una differenza significativa nella sopravvivenza tra pazienti con un buon livello di eiezione ventricolare sinistra e quelli con una frazione scarsa della stessa a partire da 3 anni dall’intervento subito.

fit.surv_lv <- prodlim(Hist(fuyrs, status) ~ lv,
     data = heart_valve
)
plot(fit.surv_lv, legend = FALSE)

Non vi è tuttavia una differenza di genere significativa su tutto il periodo di follow up.

fit.surv_lv <- prodlim(Hist(fuyrs, status) ~ sex,
     data = heart_valve
)
plot(fit.surv_lv, legend = FALSE)

La frazione di eiezione ventricolare sinistra pre-operatoria sembra presentare un impatto minore sulla sopravvivenza dei pazienti, anche se le code sembrano avere delle differenze significative in termini di probabilità di sopravvivenza. Tuttavia, va ricordato che il basso numero di pazienti a rischio, nelle code, rende le analisi non attendibili.

fit.surv_lv <- prodlim(Hist(fuyrs, status) ~ sten.reg.mix,
     data = heart_valve
)
plot(fit.surv_lv, legend = FALSE)

Suddividendo ora i soggetti in analisi in classi di età è evidente come la sopravvivenza degli over 65 sia significativamente inferiore a quella nei soggetti più giovani. Non si evidenziano invece differenze significative per i soggetti under 50 e per quelli tra i 50 e i 65 anni.

fit.surv_age <- prodlim(Hist(fuyrs, status) ~ age_lev,
     data = heart_valve
)
plot(fit.surv_age, legend = FALSE)

Soggetti con un livello di creatina elevato (\(\geq 120\)) presentano una probabilità di sopravvivenza inferiore nel lungo periodo, ovvero dopo circa 6 anni di follow up. Tuttavia la differenza non sembra essere sostanziale e la scarsità di soggetti seguiti dopo 8 anni dall’intervento alla valvola non consente di giungere a delle conclusioni.

fit.surv_creat <- prodlim(Hist(fuyrs, status) ~ creat_lev,
     data = heart_valve
)
plot(fit.surv_creat, legend = FALSE)

Anche il logaritmo dell’indice di massa ventricolare sinistra sembra presentare una relazione significativa con la sopravvivenza dei pazienti. In particolare, al crescere di tale valore si ha una probabilità di sopravvivenza decrescente.

fit.surv_creat <- prodlim(Hist(fuyrs, status) ~ log.lvmi_lev,
     data = heart_valve
)
plot(fit.surv_creat, legend = FALSE)

Infine notiamo che non è presente una forte correlazione tra le variabili continue, il che scoraggerebbe l’inserimento congiunto delle stesse in un modello predittivo. Si procede pertanto costruendo i modelli di Cox per tutte le covariate.

nums <- unlist(lapply(heart_valve, is.numeric), use.names = FALSE)
corrplot(cor(heart_valve[, nums]),
     method = "color", order = "alphabet"
)

3 Analisi univariate

Dal modello di Cox costruito per valutare l’associazione tra log.lvmi e la sopravvivenza, emerge che è presente un coefficiente \(\beta = 0.91 > 0\) e di conseguenza un aumento dell’indice di massa ventricolare sinistra peggiora l’hazard ratio, in particolare del \(48\%\). La differenza è significativa al \(5\%\) con un \(p\)-value di \(0.001\).

model1 <- coxph(
     formula = Surv(fuyrs, status) ~ log.lvmi,
     data = heart_valve
)
summary(model1)

## Call:
## coxph(formula = Surv(fuyrs, status) ~ log.lvmi, data = heart_valve)
## 
##   n= 255, number of events= 54 
## 
##            coef exp(coef) se(coef)     z Pr(>|z|)   
## log.lvmi 0.9090    2.4818   0.3514 2.587  0.00968 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##          exp(coef) exp(-coef) lower .95 upper .95
## log.lvmi     2.482     0.4029     1.246     4.942
## 
## Concordance= 0.585  (se = 0.041 )
## Likelihood ratio test= 6.64  on 1 df,   p=0.01
## Wald test            = 6.69  on 1 df,   p=0.01
## Score (logrank) test = 6.6  on 1 df,   p=0.01

Nel seguente modello si vuole testare l’associazione tra il genere e la sopravvivenza. Si nota che è presente un coefficiente \(\beta = 0.21 > 0\) e di conseguenza la variabile sex sembra aumentare il rischio di morte, nello specifico di \(1.23\). Un valore positivo, infatti, si traduce in un maggior rischio, e quindi in una prognosi peggiore, per i soggetti con valori più alti di quella variabile. Tuttavia il coefficiente non è significativo al \(5\%\).

model2 <- coxph(
     formula = Surv(fuyrs, status) ~ sex,
     data = heart_valve
)
summary(model2)

## Call:
## coxph(formula = Surv(fuyrs, status) ~ sex, data = heart_valve)
## 
##   n= 255, number of events= 54 
## 
##        coef exp(coef) se(coef)   z Pr(>|z|)
## sex1 0.2055    1.2281   0.2936 0.7    0.484
## 
##      exp(coef) exp(-coef) lower .95 upper .95
## sex1     1.228     0.8143    0.6907     2.184
## 
## Concordance= 0.526  (se = 0.035 )
## Likelihood ratio test= 0.48  on 1 df,   p=0.5
## Wald test            = 0.49  on 1 df,   p=0.5
## Score (logrank) test = 0.49  on 1 df,   p=0.5

Dall’analisi univariata dell’associazione tra la sopravvivenza e la covariata age, si osserva che il coefficiente \(\beta = 0.10\) è positivo, comportando un aumento del rischio di evento e portando ad un hazard ratio di \(11\%\). La differenza tuttavia non è significativa.

model3 <- coxph(
     formula = Surv(fuyrs, status) ~ age,
     data = heart_valve
)
summary(model3)

## Call:
## coxph(formula = Surv(fuyrs, status) ~ age, data = heart_valve)
## 
##   n= 255, number of events= 54 
## 
##        coef exp(coef) se(coef)    z Pr(>|z|)    
## age 0.10227   1.10768  0.01716 5.96 2.52e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##     exp(coef) exp(-coef) lower .95 upper .95
## age     1.108     0.9028     1.071     1.146
## 
## Concordance= 0.723  (se = 0.038 )
## Likelihood ratio test= 48.53  on 1 df,   p=3e-12
## Wald test            = 35.53  on 1 df,   p=3e-09
## Score (logrank) test = 36.68  on 1 df,   p=1e-09

Nel presente modello di Cox, implementato per testare l’associazione tra con.cabg e la sopravvivenza, si nota che è presente un coefficiente \(\beta = 0.96\) e di conseguenza la covariata peggiora l’azzardo con un hazard ratio di 2.62. Se dunque un soggetto selezionato in modo randomico ha presenza concomitante di bypass coronarico, questo ha velocità di evento maggiore di circa il \(62\%\) rispetto a chi non ce l’ha. La differenza è infatti significativa.

model4 <- coxph(
     formula = Surv(fuyrs, status) ~ con.cabg,
     heart_valve
)
summary(model4)

## Call:
## coxph(formula = Surv(fuyrs, status) ~ con.cabg, data = heart_valve)
## 
##   n= 255, number of events= 54 
## 
##             coef exp(coef) se(coef)     z Pr(>|z|)    
## con.cabg1 0.9646    2.6236   0.2731 3.532 0.000413 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##           exp(coef) exp(-coef) lower .95 upper .95
## con.cabg1     2.624     0.3812     1.536     4.481
## 
## Concordance= 0.602  (se = 0.036 )
## Likelihood ratio test= 12.03  on 1 df,   p=5e-04
## Wald test            = 12.47  on 1 df,   p=4e-04
## Score (logrank) test = 13.45  on 1 df,   p=2e-04

Con la seguente analisi si nota che è presente un coefficiente \(\beta = 0.01\), e di conseguenza la covariata peggiora l’azzardo con un hazard ratio di \(1.01\). La differenza è significativa al \(5\%\).

model5 <- coxph(
     formula = Surv(fuyrs, status) ~ creat,
     data = heart_valve
)
summary(model5)

## Call:
## coxph(formula = Surv(fuyrs, status) ~ creat, data = heart_valve)
## 
##   n= 255, number of events= 54 
## 
##           coef exp(coef) se(coef)    z Pr(>|z|)   
## creat 0.011539  1.011606 0.003821 3.02  0.00253 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##       exp(coef) exp(-coef) lower .95 upper .95
## creat     1.012     0.9885     1.004     1.019
## 
## Concordance= 0.613  (se = 0.041 )
## Likelihood ratio test= 7.7  on 1 df,   p=0.006
## Wald test            = 9.12  on 1 df,   p=0.003
## Score (logrank) test = 9.29  on 1 df,   p=0.002

Dall’analisi univariata dell’associazione tra sopravvivenza e la covariata lv, si osserva un coefficiente \(\beta = 0.49\). Un qualunque soggetto con frazione di eiezione ventricolare sinistra pre-operatoria moderata (lv=2) ha una velocità di evento maggiore del \(64\%\) rispetto ha chi ce l’ha buona (lv=1). La differenza è significativa, il valore \(1\) è incluso nell’intervallo di confidenza.

model6 <- coxph(
     formula = Surv(fuyrs, status) ~ factor(lv),
     data = heart_valve
)
summary(model6)

## Call:
## coxph(formula = Surv(fuyrs, status) ~ factor(lv), data = heart_valve)
## 
##   n= 255, number of events= 54 
## 
##               coef exp(coef) se(coef)     z Pr(>|z|)   
## factor(lv)2 0.4991    1.6473   0.2972 1.679   0.0931 . 
## factor(lv)3 1.0953    2.9901   0.4090 2.678   0.0074 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##             exp(coef) exp(-coef) lower .95 upper .95
## factor(lv)2     1.647     0.6071     0.920     2.950
## factor(lv)3     2.990     0.3344     1.341     6.665
## 
## Concordance= 0.597  (se = 0.038 )
## Likelihood ratio test= 7.13  on 2 df,   p=0.03
## Wald test            = 7.85  on 2 df,   p=0.02
## Score (logrank) test = 8.38  on 2 df,   p=0.02

Nel presente modello di Cox implementato per testare l’associazione tra sopravvivenza e l’emodinamica della valvola aortica, si nota un \(\beta = -0.86\). Concludiamo affermando che la covariata diminuisce la velocità di evento. Se un soggetto selezionato randomicamente ha come emodinamica il rigurgito ( sten.reg.mix=2), allora la sua velocità di evento è ridotta del \(37\%\) rispetto a chi ha emodinamica di stenosi. Se invece tale soggetto ha emodinamica mista (sten.reg.mix=3), ha una velocità di evento ridotta del \(32\%\) rispetto a chi ce l’ha mista. L’unica differenza significativa è quella tra la categoria 1 e 2, con un \(p\)-value di \(0.049\).

model7 <- coxph(
     formula = Surv(fuyrs, status) ~ factor(sten.reg.mix),
     data = heart_valve
)
summary(model7)

## Call:
## coxph(formula = Surv(fuyrs, status) ~ factor(sten.reg.mix), data = heart_valve)
## 
##   n= 255, number of events= 54 
## 
##                          coef exp(coef) se(coef)      z Pr(>|z|)  
## factor(sten.reg.mix)2 -0.8608    0.4228   0.4378 -1.966   0.0493 *
## factor(sten.reg.mix)3 -0.8453    0.4294   0.4371 -1.934   0.0531 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                       exp(coef) exp(-coef) lower .95 upper .95
## factor(sten.reg.mix)2    0.4228      2.365    0.1793    0.9972
## factor(sten.reg.mix)3    0.4294      2.329    0.1823    1.0114
## 
## Concordance= 0.585  (se = 0.034 )
## Likelihood ratio test= 7.75  on 2 df,   p=0.02
## Wald test            = 6.75  on 2 df,   p=0.03
## Score (logrank) test = 7.17  on 2 df,   p=0.03

4 Modelli multivariati di Cox

4.1 Modello di Cox base

Lo sviluppo del modello predittivo base prevede l’inserimento dei seguenti regressori: sex, age, con.cabg, creat, lv, sten.reg.mix

Le variabili significative, date dal modello di Cox base sono: - L’età - La presenza concomitante di bypass coronarico - La frazione di eiezione ventricolare sinistra pre-operatoria I loro coefficienti sono tutti positivi per cui il rischio di morte è maggiore, e quindi la prognosi peggiore, per i soggetti con valori più alti di quelle variabili. Un aumento di età, un aumento della presenza concomitante di bypass coronarico e di una scarsa frazione di eiezione ventricolare sinistra pre-operatoria accrescono la velocità di evento. Per esempio, al netto delle altre variabili, abbiamo un aumento del \(10\%\) di hazard ad ogni anno in più di età.

model_base <- coxph(
     formula = Surv(fuyrs, status) ~ sex + age + con.cabg + creat + factor(lv) + factor(sten.reg.mix),
     data = heart_valve,
     x = T
)
summary(model_base)

## Call:
## coxph(formula = Surv(fuyrs, status) ~ sex + age + con.cabg + 
##     creat + factor(lv) + factor(sten.reg.mix), data = heart_valve, 
##     x = T)
## 
##   n= 255, number of events= 54 
## 
##                            coef exp(coef)  se(coef)      z Pr(>|z|)    
## sex1                   0.275329  1.316964  0.321802  0.856 0.392228    
## age                    0.098936  1.103996  0.019485  5.078 3.82e-07 ***
## con.cabg1              0.650285  1.916086  0.295821  2.198 0.027932 *  
## creat                  0.006730  1.006753  0.004461  1.509 0.131411    
## factor(lv)2            0.600436  1.822913  0.306889  1.957 0.050403 .  
## factor(lv)3            1.645095  5.181502  0.446836  3.682 0.000232 ***
## factor(sten.reg.mix)2 -0.619062  0.538449  0.468112 -1.322 0.186013    
## factor(sten.reg.mix)3 -0.812285  0.443843  0.448413 -1.811 0.070069 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                       exp(coef) exp(-coef) lower .95 upper .95
## sex1                     1.3170     0.7593    0.7009     2.475
## age                      1.1040     0.9058    1.0626     1.147
## con.cabg1                1.9161     0.5219    1.0730     3.422
## creat                    1.0068     0.9933    0.9980     1.016
## factor(lv)2              1.8229     0.5486    0.9989     3.327
## factor(lv)3              5.1815     0.1930    2.1583    12.439
## factor(sten.reg.mix)2    0.5384     1.8572    0.2151     1.348
## factor(sten.reg.mix)3    0.4438     2.2530    0.1843     1.069
## 
## Concordance= 0.784  (se = 0.034 )
## Likelihood ratio test= 70.28  on 8 df,   p=4e-12
## Wald test            = 48.32  on 8 df,   p=9e-08
## Score (logrank) test = 56.78  on 8 df,   p=2e-09

4.2 Valutazione dell’ipotesi PH

Dalla verifica dell’assunzione di linearità della variabile Age sembra che l’aumento della velocità di evento con l’accrescere dell’età dei pazienti sia costante nel tempo. Infatti, attraverso la rappresentazione grafica dei Martingale Residuals, non si osserva un trend persistente e si può quindi supporre che l’effetto sull’hazard sia costante all’aumentare dell’età. Abbiamo ha conferma di questo fatto anche dalla visualizzazione del modello di Cox implementato con una b-spline per la variabile age: anche in questo caso si osserva che l’età ha un effetto pressoche lineare sui log-hazard e che la curva ha un andamento complessivamente costante man mano che aumenta la densità della variabile age.

par(mfrow = c(2, 1), mar = c(4, 4, 2, 2))

mar.res <- resid(model_base, type = "martingale")

plot(heart_valve$age, mar.res,
     xlab = "age", ylab = "Martingale Residuals",
     main = "Check functional form of age"
)

lines(lowess(heart_valve$age, mar.res), col = "red")
model.age.bs <- coxph(Surv(fuyrs, status) ~ bs(age, 4),
     data = heart_valve
)

par(mar = c(4, 4, 1, 1))

plotHR(model.age.bs,
     term = "age", plot.bty = "o", ylog = T, xlim = c(30, 100),
     rug = "density", xlab = "age", polygon_ci = T
)

Per quanto riguarda il controllo sulla linearità della variabile creat, sembra non ci sia un andamento complessivamente lineare dei residui. Ciò potrebbe essere dovuto al fatto che nel modello di Cox, questa non sia una variabile significativa. Tuttavia, dalla rappresentazione grafica Martingale Residuals, non si osserva un netto trend dei residui quindi si può supporre che l’effetto sull’hazard sia tendenzialmente costante all’aumentare della creatina. Osservando meglio la forma funzionale della variabile continua sul grafico del modello di Cox implementato con una b-spline per la variabile creat, si osserva un andamento lineare della curva per valori alti di densità della variabile.

par(mfrow = c(2, 1), mar = c(4, 4, 2, 2))

mar.res <- resid(model_base, type = "martingale")
plot(heart_valve$creat, mar.res,
     xlab = "creat", ylab = "Martingale Residuals",
     main = "Check functional form of creatine"
)

lines(lowess(heart_valve$creat, mar.res), col = "green")
model.creat.bs <- coxph(Surv(fuyrs, status) ~ bs(creat, 4),
     data = heart_valve
)

par(mar = c(4, 4, 1, 1))

plotHR(model.creat.bs,
     term = "creat", plot.bty = "o", ylog = T, xlim = c(30, 100),
     rug = "density", xlab = "creat", polygon_ci = T
)

Dai grafici per i residui di Schoenfeld sembra si possa assumere PH per la variabile sex in quanto l’andamento dei residui è orizzontale, mentre necessita un’ulteriore controllo sulla linearità dei residui la variabile con.cabg. Tuttavia, osservando i \(p\)-value per sex (\(0.71\)) e per con.cabg (\(0.37\)), sembra che ciascuna delle variabili non siano significative (\(p\)-value > \(\alpha\)) e che quindi si possa assumere PH per entrambe.

par(mfrow = c(2, 1), mar = c(4, 4, 2, 2))

checkPH.sex <- cox.zph(model_base)[1]
plot(checkPH.sex, main = "Check PH assumption of sex")
points(checkPH.sex$x, checkPH.sex$y, pch = 16, col = "red")
abline(h = 0, lty = 2, col = 2)

checkPH.con_cabg <- cox.zph(model_base)[3]
plot(checkPH.con_cabg, main = "Check PH assumption of con.cabg")
points(checkPH.con_cabg$x, checkPH.con_cabg$y, pch = 16, col = "red")
abline(h = 0, lty = 2, col = 2)

Dal controllo grafico della PH assumption implementato calcolando la stima K-M nei diversi livelli della variabile categorica con.cabg, si può assumere PH per tale variabile. Nonostante i logaritmi \(\log(-\log(S_X(t)))\) e \(\log(-\log(S_0(t)))\) non si mantengano ad una distanza costante nel tempo, non sussiste un cambio di inversione di pendenza per le due classi nel tempo.

km.con.cabg <- survfit(Surv(fuyrs, status) ~ con.cabg,
     data = heart_valve
)
plot(km.con.cabg,
     col = c("black", "red"), fun = "cloglog",
     ylab = "log(-log(Survival))", xlab = "log(time)",
     main = "Check PH assumption of con.cabg"
)

Anche in questo caso la rappresentazione dei residui di Schoenfeld ci permette di assumere PH sia per lv che per sten.reg.mix: i residui hanno trend orizzontale e non dipendono dal tempo.

par(mfrow = c(2, 1), mar = c(4, 4, 2, 2))

checkPH.lv <- cox.zph(model_base)[5]
plot(checkPH.lv, main = "Check PH assumption of lv")
points(checkPH.lv$x, checkPH.lv$y, pch = 16, col = "red")
abline(h = 0, lty = 2, col = 2)

checkPH.sten <- cox.zph(model_base)[6]
plot(checkPH.sten, main = "Check PH assumption of sten reg")
points(checkPH.sten$x, checkPH.sten$y, pch = 16, col = "red")
abline(h = 0, lty = 2, col = 2)

4.3 Modello di Cox aumentato

Lo sviluppo del modello predittivo aumentato prevede l’inserimento dei seguenti regressori: sex, age, con.cabg, creat, lv, sten.reg.mix e log.lvmi. Come descritto in precedenza, il modello di Cox aumentato differisce dal modello base per il solo inserimento aggiuntivo del logaritmo naturale dell’indice di massa ventricolare sinistra. La statistica del test di Wald è superiore a quella del modello base, pertanto da un’analisi preliminare sembra esserci un miglioramento.

model_augmented <- coxph(
     formula = Surv(fuyrs, status) ~ log.lvmi + sex + age + con.cabg + creat + lv + sten.reg.mix,
     data = heart_valve, x = T
)
summary(model_augmented)

## Call:
## coxph(formula = Surv(fuyrs, status) ~ log.lvmi + sex + age + 
##     con.cabg + creat + lv + sten.reg.mix, data = heart_valve, 
##     x = T)
## 
##   n= 255, number of events= 54 
## 
##                    coef exp(coef)  se(coef)      z Pr(>|z|)    
## log.lvmi       1.219381  3.385092  0.388812  3.136 0.001712 ** 
## sex1           0.395847  1.485642  0.326592  1.212 0.225493    
## age            0.098178  1.103159  0.019689  4.987 6.15e-07 ***
## con.cabg1      0.737046  2.089753  0.300918  2.449 0.014313 *  
## creat          0.003480  1.003486  0.004706  0.739 0.459702    
## lv2            0.617825  1.854889  0.311077  1.986 0.047024 *  
## lv3            1.662428  5.272097  0.445345  3.733 0.000189 ***
## sten.reg.mix2 -0.681502  0.505857  0.474146 -1.437 0.150625    
## sten.reg.mix3 -0.935386  0.392434  0.448925 -2.084 0.037195 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##               exp(coef) exp(-coef) lower .95 upper .95
## log.lvmi         3.3851     0.2954    1.5798     7.253
## sex1             1.4856     0.6731    0.7833     2.818
## age              1.1032     0.9065    1.0614     1.147
## con.cabg1        2.0898     0.4785    1.1587     3.769
## creat            1.0035     0.9965    0.9943     1.013
## lv2              1.8549     0.5391    1.0082     3.413
## lv3              5.2721     0.1897    2.2024    12.620
## sten.reg.mix2    0.5059     1.9768    0.1997     1.281
## sten.reg.mix3    0.3924     2.5482    0.1628     0.946
## 
## Concordance= 0.796  (se = 0.032 )
## Likelihood ratio test= 79.86  on 9 df,   p=2e-13
## Wald test            = 58.92  on 9 df,   p=2e-09
## Score (logrank) test = 67.18  on 9 df,   p=5e-11

Dai test di Wald sulle singole covariate risultano significative per spiegare l’azzardo: - Il logaritmo dell’indice di massa ventricolare - L’età - La presenza di bypass coronario concomitante. Inoltre, al crescere del logaritmo di lvmi, al netto delle altre covariate, si ha un aumento dell’hazard ratio e di conseguenza un peggioramento della prognosi. Lo stesso vale per l’età: individui più anziani avranno una prognosi peggiore e una conseguente probabilità di sopravvivenza minore rispetto ad individui più giovani. 2 Inoltre, individui che presentano un bypass coronario concomitante avranno un hazard ratio superiore e una probabilità di sopravvivenza più bassa.

Si sottolinea inoltre che vi è una netto contributo alla previsione dell’azzardo dato dalla frazione di eiezione ventricolare sinistra e dall’emodinamica della valvola aortica. Nello specifico, un individuo con un livello moderato o scarso della frazione di eiezione ventricolare sinistra avrà una probabilità di sopravvivenza nel tempo inferiore rispetto ad un individuo con un buon livello, a parità di altri fattori. Viceversa, una valvola aortica mista rispetto ad una valvola stenotica ha un effetto negativo sull’azzardo e conseguentemente, a parità degli altri fattori, la prognosi per un individuo con valvola aortica mista sarà migliore.

Il genere e il livello di creatina serica non risultano invece significativi nello spiegare la variabilità dell’azzardo.

4.4 Valutazione dell’ipotesi PH

Dal seguente test sui residui di Schoenfeld, l’ipotesi di Proportional Hazard sembra verificata per tutte le variabili.

par(mfrow = c(2, 1), mar = c(4, 4, 2, 2))

checkPH.kar <- cox.zph(model_augmented)
checkPH.kar

##                chisq df     p
## log.lvmi      6.2457  1 0.012
## sex           0.0094  1 0.923
## age           1.7535  1 0.185
## con.cabg      0.6038  1 0.437
## creat         0.0467  1 0.829
## lv            0.0529  2 0.974
## sten.reg.mix  1.0247  2 0.599
## GLOBAL       12.2951  9 0.197

Per quanto riguarda log.lvmi, invece, si rifiuterebbe l’ipotesi alternativa di azzardo non proporzionale nel tempo solo per un livello di confidenza del \(99\%\). Si procede quindi con l’analisi grafica delle distribuzioni dei residui di Schoenfeld per le variabili categoriali e dei residui di Martingale per le variabili continue.

La stima dei parametri \(\mathbf{\beta}\) tempo dipendenti per le variabili genere, frazione di eiezione ventricolare sinistra ed emodinamica della valvola aortica, evidenziano un andamento approssimativamente costante nel tempo. Si ha pertanto un’ulteriore conferma circa la validità dell’assunzione di azzardo proporzionale per queste variabili. Per quanto concerne il parametro tempo dipendente del modello di Cox relativo alla presenza concomitante di bypass coronario, si osserva un andamento parabolico che porta a dubitare di un azzardo proporzionale nel tempo.

# sex
par(mfrow = c(2, 2))

plot(checkPH.kar[2], main = "Check PH assumption of sex")
points(checkPH.kar[2]$x, checkPH.kar[2]$y, pch = 16, col = "lightgray")
abline(h = 0, lty = 2, col = 2)

# con.cabg
plot(checkPH.kar[4], main = "Check PH assumption of con.cabg")
points(checkPH.kar[4]$x, checkPH.kar[4]$y, pch = 16, col = "green")
abline(h = 0, lty = 2, col = 2)

# sten.reg.mix
plot(checkPH.kar[7], main = "Check PH assumption of sten.reg.mix")
points(checkPH.kar[7]$x, checkPH.kar[7]$y, pch = 16, col = "red")
abline(h = 0, lty = 2, col = 2)

# lv
plot(checkPH.kar[6], main = "Check PH assumption of lv")
points(checkPH.kar[6]$x, checkPH.kar[6]$y, pch = 16, col = "blue")
abline(h = 0, lty = 2, col = 2)

Tuttavia, tramite analisi grafica della differenza dei logaritmi delle due curve di sopravvivenza, possiamo assumere che valga PH anche in questo caso

km.con.cabg <- survfit(Surv(fuyrs, status) ~ strata(con.cabg), data = heart_valve)
plot(km.con.cabg,
     col = c("black", "red"), fun = "cloglog",
     ylab = "log(-log(Survival))", xlab = "log(time)",
     main = "Check PH assumption of con.cabg"
)

I residui di Martingale per il logaritmo del lvmi presentano una media decrescente nel tempo. Inoltre la forma funzionale dell’hazard ratio per tale variabile non è lineare, ma parabolica. Non potendo confermare l’ipotesi di azzardo proporzionale nel tempo, si costruisce un modello con una variabile tempo-dipendente al fine di verificare se il logaritmo di lvmi ha un impatto sull’azzardo dipendente dalla componente temporale.

km.log.lvmi <- coxph(Surv(fuyrs, status) ~ log.lvmi,
     data = heart_valve
)

par(mfrow = c(2, 1), mar = c(4, 4, 2, 2))

mar.res <- resid(km.log.lvmi, type = "martingale", times = )
plot(heart_valve$log.lvmi, mar.res,
     xlab = "lvmi", ylab = "Martingale Residuals",
     main = "Check functional form of log lvmi"
)
lines(lowess(heart_valve$log.lvmi, mar.res), col = "red")

model.kar.bs <- coxph(Surv(fuyrs, status == 1) ~ bs(log.lvmi, 4),
     data = heart_valve
)
plotHR(model.kar.bs,
     term = "log.lvmi", plot.bty = "o",
     xlim = c(0, 10), xlab = "lvmi"
)

Viene infine inserita una variabile tempo dipendente per log.lvmi in modo da valutare l’assunto di azzardo proporzionale messo in dubbio precedentemente. La variabile tempo-dipendente è data dal prodotto di log.lvmi e del logaritmo di un terzo del valore del tempo: \[ g(t) = \log(t) - \log(3) \] Il test di Wald per tale variabile non risulta però significativo, suggerendo un azzardo proporzionale nel tempo, come da ipotesi del modello di Cox.

cut.points <- unique(heart_valve$fuyrs[heart_valve$status == 1])
heart_valve$entry <- 0
heart_valve2 <- survSplit(
     data = heart_valve, cut = cut.points,
     end = "fuyrs", start = "entry", event = "status"
)
heart_valve3 <- heart_valve2[order(heart_valve2$paz.id), ]
heart_valve3$td_lvmi <- heart_valve3$log.lvmi * log(heart_valve3$fuyrs / 3)
modello_td <- coxph(Surv(entry, fuyrs, status) ~ log.lvmi + sex + age + con.cabg + creat + lv + sten.reg.mix + td_lvmi,
     data = heart_valve3
)
summary(modello_td)

## Call:
## coxph(formula = Surv(entry, fuyrs, status) ~ log.lvmi + sex + 
##     age + con.cabg + creat + lv + sten.reg.mix + td_lvmi, data = heart_valve3)
## 
##   n= 8524, number of events= 54 
## 
##                    coef exp(coef)  se(coef)      z Pr(>|z|)    
## log.lvmi       1.256452  3.512936  0.394038  3.189 0.001429 ** 
## sex1           0.420603  1.522879  0.327790  1.283 0.199441    
## age            0.098542  1.103560  0.019693  5.004 5.62e-07 ***
## con.cabg1      0.737683  2.091085  0.301040  2.450 0.014268 *  
## creat          0.003754  1.003761  0.004726  0.794 0.426966    
## lv2            0.632255  1.881849  0.311502  2.030 0.042387 *  
## lv3            1.678120  5.355478  0.445372  3.768 0.000165 ***
## sten.reg.mix2 -0.686334  0.503418  0.474050 -1.448 0.147670    
## sten.reg.mix3 -0.959949  0.382912  0.450022 -2.133 0.032915 *  
## td_lvmi        0.273060  1.313979  0.366670  0.745 0.456451    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##               exp(coef) exp(-coef) lower .95 upper .95
## log.lvmi         3.5129     0.2847    1.6228     7.605
## sex1             1.5229     0.6567    0.8010     2.895
## age              1.1036     0.9062    1.0618     1.147
## con.cabg1        2.0911     0.4782    1.1591     3.772
## creat            1.0038     0.9963    0.9945     1.013
## lv2              1.8818     0.5314    1.0220     3.465
## lv3              5.3555     0.1867    2.2372    12.820
## sten.reg.mix2    0.5034     1.9864    0.1988     1.275
## sten.reg.mix3    0.3829     2.6116    0.1585     0.925
## td_lvmi          1.3140     0.7610    0.6404     2.696
## 
## Concordance= 0.795  (se = 0.033 )
## Likelihood ratio test= 80.45  on 10 df,   p=4e-13
## Wald test            = 58.95  on 10 df,   p=6e-09
## Score (logrank) test = 67.31  on 10 df,   p=1e-10

I residui di Martingale per l’età si distribuiscono in modo approssimativamente casuale attorno alla media, con una leggera differenza per individui molto anziani. La distribuzione dell’Hazard Ratio ottenuta con il metodo di Karnofsky evidenzia linearità, tenendo anche in considerazione la distribuzione d’età del campione di riferimento. L’ipotesi di azzardo proporzionale per l’età è dunque confermata.

km.age <- coxph(Surv(fuyrs, status) ~ age,
     data = heart_valve
)

par(mfrow = c(2, 1), mar = c(4, 4, 2, 2))

mar.res_age <- resid(km.age, type = "martingale", times = )
plot(heart_valve$age, mar.res_age,
     xlab = "age", ylab = "Martingale Residuals",
     main = "Check functional form of age"
)
lines(lowess(heart_valve$age, mar.res_age), col = "red")

model.age.bs <- coxph(Surv(fuyrs, status == 1) ~ bs(age, 4),
     data = heart_valve
)
plotHR(model.age.bs,
     term = "age", plot.bty = "o",
     xlim = c(0, 400), xlab = "age"
)

Anche i residui di Martingale per la creatina si distribuiscono in modo approssimativamente casuale attorno alla media. Inoltre l’andamento dell’Hazard Ratio nell’intervallo di livello di creatina ad alta concentrazione è indicativamente lineare. Si ha pertanto evidenza di azzardo proporzionale anche per il livello di creatina serica pre-operatoria.

km.creat <- coxph(Surv(fuyrs, status) ~ creat,
     data = heart_valve
)

par(mfrow = c(2, 1), mar = c(4, 4, 2, 2))

mar.res_creat <- resid(km.creat, type = "martingale", times = )
plot(heart_valve$creat, mar.res_creat,
     xlab = "creat", ylab = "Martingale Residuals",
     main = "Check functional form of creatine"
)
lines(lowess(heart_valve$creat, mar.res_creat), col = "red")

model.creat.bs <- coxph(Surv(fuyrs, status == 1) ~ bs(creat, 4),
     data = heart_valve
)
plotHR(model.creat.bs,
     term = "creat", plot.bty = "o",
     xlim = c(0, 400), xlab = "creat"
)

5 Validazione e confronto

Per iniziare, istanziamo il numero di anni di cui siamo interessati a prevedere il rischio. In questo caso è 5.

N <- 5

Calcoliamo poi la sopravvivenza predetta dai due modelli, e stimiamo il rischio di morte entro 5 anni per ogni paziente.

# Sopravvivenza di ogni paziente, per ogni modello
fit_base <- survfit(model_base, newdata = heart_valve)
fit_augmented <- survfit(model_augmented, newdata = heart_valve)

# Rischio di morte entro 5 anni
heart_valve$risk_of_death_base <- 1 - as.numeric(summary(fit_base, times = N)$surv)
heart_valve$risk_of_death_augmented <- 1 - as.numeric(summary(fit_augmented, times = N)$surv)

Ora, mancano da calcolare gli score necessari per la costruzione di indici grafici di valutazione dei modelli. In particolare, siamo interessati alla loro calibrazione e discriminazione.
Per quanto riguarda la calibrazione, testiamo in particolare la Weak Calibration, sfruttando il Calibration Plot e il Brier Score. Tale indice numerico è successivamente confrontato con lo stesso indice calcolato nell’ipotesi di Strong Calibration.
Per quanto riguarda invece la discriminazione, osserviamo il comportamento delle due curve ROC. Da esse si può poi estrarre l’indice di area sotto la curva, ovvero l’AUC, e valutarne la differenza. Infine, esploriamo anche come le AUC si comportano nel tempo per ognuno dei due modelli.

Sfruttiamo quindi la libreria riskRegression per calcolare gli score necessari alle analisi. Calcoliamo questi score per ogni tempo, da 1 a 5. Non contiamo i tempi successivi, per confinare la nostra domanda di ricerca ad un orizzonte temporale limitato, ma anche per la poca affidabilità dovuta alla numerosità campionaria.
Si noti che si non si usa un procedimento di Cross-Validation, incompatibile con il calcolo della ROC curve nel pacchetto riskRegression.

score <- Score(list("base" = model_base, "augmented" = model_augmented),
     formula = Surv(fuyrs, status) ~ 1,
     data = heart_valve,
     times = 1:N, plots = c("cal", "roc", "auc"), summary = "risks"
)

5.1 Calibrazione

Dal Calibration Plot notiamo comportamenti simili tra i due modelli, in termini di rischi stimati. Benchè correlati tra loro, entrambi sovrastimano il rischio nelle classi di pazienti selezionate da alte soglie di rischio. In generale, l’atteggiamento dei modelli è quindi di prudenza: d’altra parte avremmo uno svantaggio per i pazienti se i modelli sottostimassero il rischio.
Da un’analisi più approfondita si nota che precedentemente al valore di soglia del \(35\%\) il modello aumentato sottostima il rischio, mentre dopo tale soglia lo sovrastima.

# Calibration plot at 5 years
par(mfrow = c(1, 1), mar = c(4, 4, 2, 2))

plotCalibration(score,
     times = N, type = "p",
     auc.in.legend = F, brier.in.legend = F,
     cens.method = "local",
     col = c("#08519c", "#86bff9")
)

Notiamo che anche i Brier Score sono coerenti tra loro, con un valore attorno a \(0.35\). Tralasciando il loro valore assoluto ma confrontandoli con i Brier Score calcolati nell’ipotesi di Strong Calibration, osserviamo che sono molto simili tra loro, seppur non approssimabili a zero.
Il modello che in generale ha la migliore calibrazione è il modello base, con previsioni più vicine alla diagonale del quadrante del Calibration Plot.

# Indicatore di evento entro 5 anni
heart_valve$event.Ny <- ifelse(heart_valve$fuyrs <= N, 1, 0)

# Brier Score
BS_base <- mean((heart_valve$event.Ny - heart_valve$risk_of_death_base)^2)
BS_augmented <- mean((heart_valve$event.Ny - heart_valve$risk_of_death_augmented)^2)
print(c("BS_base: ", round(BS_base, digits = 3)))

## [1] "BS_base: " "0.348"

print(c("BS_augmented: ", round(BS_augmented, digits = 3)))

## [1] "BS_augmented: " "0.355"

# Brier Score sotto l'ipotesi di Strong Calibration
BS_base_sc <- mean(heart_valve$risk_of_death_base * (1 - heart_valve$risk_of_death_base))
BS_augmented_sc <- mean(heart_valve$risk_of_death_augmented * (1 - heart_valve$risk_of_death_augmented))

# Si comparano i Brier Score osservati con quelli attesi sotto Strong Calibration
print(c("BS_base - BS_base_sc: ", round(BS_base - BS_base_sc, digits = 3)))

## [1] "BS_base - BS_base_sc: " "0.245"

print(c("BS_augmented - BS_augmented_sc: ", round(BS_augmented - BS_augmented_sc, digits = 3)))

## [1] "BS_augmented - BS_augmented_sc: " "0.26"

5.2 Discriminazione

Proseguiamo le analisi tramite curve ROC e Youden Index.
La prudenza del modello aumentato notata nel Calibration Plot si traduce qui in una curva ROC superiore a quella del modello base solo per soglie di rischio meno stringenti, e quindi per una frazione di popolazione più ampia. Il modello base, al contrario, ha una miglior curva per soglie più stringenti.
Le due AUC, seppur vicine, sottolineano comunque un leggero vantaggio per il modello aumentato.

# ROC curve a 5 anni
plotROC(score,
     times = N, xlab = "FPR", ylab = "TPR",
     col = c("#08519c", "#86bff9")
)

Calcoliamo ora le soglie ottimali secondo lo Youden Index. Notiamo che la soglia ottima del modello aumentato è, come ci aspettiamo, più ampia: sceglie di classificare come a rischio i pazienti con un rischio superiore all’\(11\%\), in confronto al \(30\%\) del modello base.

# Dati delle curve ROC di ogni modello
roc <- score$ROC$plotframe |> data.frame()
roc <- roc[roc$times == N, ]
roc_base <- roc[roc$model == "base", ]
roc_augmented <- roc[roc$model == "augmented", ]

# Youden indices
Youden_base <- roc_base$TPR - roc_base$FPR
Youden_augmented <- roc_augmented$TPR - roc_augmented$FPR

# Soglie ottimali
opt_cutoff_base <- roc_base$risk[Youden_base == max(Youden_base)]
opt_cutoff_augmented <- roc_augmented$risk[Youden_augmented == max(Youden_augmented)]

print(c("Optimal cut off base model ", round(opt_cutoff_base, digits = 3)))

## [1] "Optimal cut off base model " "0.305"

print(c("Optimal cut off augmented ", round(opt_cutoff_augmented, digits = 3)))

## [1] "Optimal cut off augmented " "0.107"

# Plot delle soglie sopra le curve ROC
plotROC(score,
     times = N, xlab = "FPR", ylab = "TPR",
     col = c("#08519c", "#86bff9")
)

points(roc_base$FPR[Youden_base == max(Youden_base)],
     roc_base$TPR[Youden_base == max(Youden_base)],
     pch = 22, cex = 1.7, bg = "#08519c"
)

points(roc_augmented$FPR[Youden_augmented == max(Youden_augmented)],
     roc_augmented$TPR[Youden_augmented == max(Youden_augmented)],
     pch = 22, cex = 1.7, bg = "#86bff9"
)

La curva AUC relativa al secondo modello sovrasta, di poco, quella del primo modello: ne segue dunque una maggior capacità discriminante all’aumentare dell’orizzonte temporale.

# AUC
plotAUC(score,
     which = "score",
     col = c("#08519c", "#86bff9")
)

Proviamo ora ad interpretare il Net Benefit Plot.
E’ proprio intorno ad un \(35\%\) di soglia ottimale, che troviamo il beneficio maggiore, a favore del modello base. Intorno a questa soglia abbiamo un plateau delle curve di beneficio, con un vantaggio per il modello base. Il modello aumentato appare migliore a basse soglie, coerentemente con le precedenti analisi.

# Net benefit al tempo 5
dca(Surv(fuyrs, status) ~ risk_of_death_base + risk_of_death_augmented,
     data = heart_valve, time = N
)

## Printing with `plot(x, type = 'net_benefit', smooth = FALSE, show_ggplot_code = FALSE)`

In conclusione, la scelta del modello clinico da adottare può dipendere da vari fattori. Per una maggiore prudenza, che si traduce in una sovrastima del rischio e soglie ottimali più basse, si preferisce il modello aumentato. Adottando questo modello, si mette in secondo piano il costo di miclassificare pazienti non a rischio come a rischio.
Se, al contrario, si mette al primo piano il costo di misclassificazione ed il valore di ottenere un modello calibrato e preciso, si scelga il modello base.

Seppur non riportati nel presente lavoro, gli stessi grafici calcolati con un orizzonte temporale più ampio (precisamente di 8 anni), avvallano la medesima tesi.

6 Previsione del rischio

Si vuole prevedere il rischio di evento ad un time-point fisso, per tre soggetti tipo. Per l’individuazione di 3 soggetti tipo, si fa affidamento ad una tecnica di clustering. L’idea di base dell’algoritmo K-means si traduce, nel caso di variabili miste, nell’algoritmo K-medoids. Si cerca dunque di clusterizzare il sample di pazienti in 3 gruppi, e di estrarre successivamente le osservazioni centrali.
La metrica di distanza usata è quella di Gower, come suggerito dalla letteratura.

# Data matrix per il clustering, escludendo le variabili: id paziente, tempo di evento e status
data_matrix <- heart_valve[, c(2, 5, 6, 7, 8, 9, 10)]

# Matrice di distanza (gower distance)
gower_distance_matrix <- as.matrix(daisy(data_matrix, metric = "gower"))

# K-medoids clustering, con 3 centri
n_centers <- 3
clusters <- pam(gower_distance_matrix, n_centers,
     metric = "euclidean",
     medoids = "random",
     nstart = 10
)

# Estraiamo i centri dall'output del clustering
subjects_id <- as.vector(clusters$id.med)
subjects <- data_matrix[subjects_id, ]

I seguenti risultati rappresentano i rischi dei 3 soggetti a 5 anni. Il modello base ottiene rischi nettamente maggiori.

# Stimiamo il rischio di morte dei 3 soggetti a 5 anni
for (i in 1:n_centers) {
     fit1 <- survfit(model_base, newdata = subjects[i, ])$surv
     fit2 <- survfit(model_augmented, newdata = subjects[i, ])$surv

     risk_N_base <- 1 - fit1[249 / max(heart_valve$fuyrs) * N]
     risk_N_augmented <- 1 - fit2[249 / max(heart_valve$fuyrs) * N]

     print(c("ID:", subjects_id[i]))
     print(c("Base model risk: ", round(risk_N_base, digits = 3)))
     print(c("Augmented model risk: ", round(risk_N_augmented, digits = 3)))
}

## [1] "ID:" "216"
## [1] "Base model risk: " "0.217"            
## [1] "Augmented model risk: " "0.179"                 
## [1] "ID:" "218"
## [1] "Base model risk: " "0.21"             
## [1] "Augmented model risk: " "0.113"                 
## [1] "ID:" "9"  
## [1] "Base model risk: " "0.007"            
## [1] "Augmented model risk: " "0.003"

Provando anche a visualizzare l’evoluzione del rischio nei tre soggetti, per ognuno dei due modelli, e confrontandolo con la curva di sopravvivenza del dataset, possiamo notare come le previsioni del modello aumentato diano una sopravvivenza maggiore. Questo è in linea con gli indici di rischio a 5 anni, trovati in precedenza.

par(mfrow = c(2, 1), mar = c(4, 4, 2, 2))

# Curva di sopravvivenza osservata
general_surv <- survfit(Surv(fuyrs, status) ~ 1,
     data = heart_valve
)$surv

# Modello base: previsione per i 3 soggetti
plot(general_surv,
     type = "s", col = 1, ylim = c(0.6, 1), xlim = c(0, 185), lty = 2, lwd = 2,
     xlab = "time", ylab = "survival probability"
)
for (i in 1:n_centers) {
     fit <- survfit(model_base, newdata = subjects[i, ])$surv
     lines(fit, type = "s", lwd = 1.5, col = i + 4)
}

# Modello aumentato: previsione per i 3 soggetti
plot(general_surv,
     type = "s", col = 1, ylim = c(0.6, 1), xlim = c(0, 185), lty = 2, lwd = 2,
     xlab = "time", ylab = "survival probability"
)
for (i in 1:n_centers) {
     fit <- survfit(model_augmented, newdata = subjects[i, ])$surv
     lines(fit, type = "s", lwd = 1.5, col = i + 4)
}